/* Programa que realiza a multiplicação de uma matriz por um vetor */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#include <string.h>
#include <unistd.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#include <gsl/gsl_rng.h>

#include <time.h>

// Modifiquei




void imprimeMatriz(FILE* fp, gsl_matrix* P, int dim);
gsl_matrix* calcProduto(gsl_matrix* A, gsl_matrix* B, int dim);
void initMatrizes(gsl_rng * r, gsl_matrix* A, int dim);

int main(void)
{
	int N = 10000; 
	int dim = 2;
	int i;
	FILE *fp;
	const gsl_rng_type * T;
	gsl_rng * r;
	
	gsl_rng_env_setup();
	
	T = gsl_rng_default;
	r = gsl_rng_alloc(T);
	
	gsl_rng_set(r, (int)time(NULL));
	
	gsl_matrix* A[N]; /* As matrizes aleatórias são guardadas em A */ 
	gsl_matrix* P[N]; /* Os produtos parciais são guardados em P */ 
	
	fp = fopen("resultado-mult", "w");
	
	for ( i = 0 ; i < N ; i++ )
	{
		A[i] = gsl_matrix_alloc(dim, dim);
		P[i] = gsl_matrix_alloc(dim, dim);
	}
	
	for ( i = 0 ; i < N ;i++ )
	{
		initMatrizes(r, A[i], dim);
		imprimeMatriz(fp, A[i], dim);
	}
	
	
	
	/* Primeiro Produto */
	
	P[0] = calcProduto(A[0], A[1], dim);
	imprimeMatriz(fp, P[0], dim);
	
	// Produtos parciais
	for ( i = 1 ; i < (N-1) ; i++ )
	{
		P[i] = calcProduto(P[i-1], A[i+1], dim);
		imprimeMatriz(fp, P[i], dim);
	}

	//imprimeMatriz(fp, P[N-2], dim);
	
	fclose(fp);
	return 0;

}

void imprimeMatriz(FILE* fp, gsl_matrix* P, int dim)
{
	int i,j;
	
	for ( i = 0 ; i < dim ; i++ )
	{
		for ( j = 0 ; j < dim ; j++ )
		{
			fprintf(fp, "%14.6lf\t", gsl_matrix_get(P, i, j));
		}
		fprintf(fp, "\n");
	}
	fprintf(fp, "\n");
	return;
}

gsl_matrix* calcProduto(gsl_matrix* A, gsl_matrix* B, int dim)
{
	gsl_matrix* R;
	double tmp;
	int i,j,k;
	
	R = gsl_matrix_alloc(dim, dim);
	
	gsl_matrix_set_zero(R);
	
	for ( i = 0 ; i < dim ; i++ )
	{
		for ( j = 0 ; j < dim ; j++ )
		{
			tmp = 0;
			for ( k = 0 ; k < dim ; k++ )
			{
				tmp = tmp + gsl_matrix_get(A, i, k)*gsl_matrix_get(B, k, j);
			}
			gsl_matrix_set(R, i, j, tmp);
		}
	}
	
	return R;
}

void initMatrizes(gsl_rng * r, gsl_matrix* A, int dim)
{
	int i,j;
	double num;
	static double total = 0;
	
	for ( i = 0 ; i < dim ; i++ )
	{
		for ( j = 0 ; j < dim ; j++ )
		{
			num = gsl_rng_uniform_pos(r);
			gsl_matrix_set(A, i, j, num);
			total = total + num;
		}
	}
	
	printf("Total = %f\n", total/40000);
	
	return;
}

/* gcc -lm -lgsl -O3 -lgslcblas mult_matrizes_seq.c -o mult_matrizes_seq.exe */